function [L,U] = doolittle_decomposition(X)
L = zeros(size(X));
U = zeros(size(X));
for j = 1:length(X)
    U(1,j) = X(1,j);
end
for j = 1:length(X)
    L(j,1) = X(j,1)/X(1,1);
    L(j,j) = 1;
end
for i = 2:length(X)
    for j = 2:length(X)
        if i == j
            s = 0;
            for k = 1:i-1
                s = s + L(i,k)*U(k,i);
            end
            U(i,i) = X(i,i) - s;
        elseif i > j
            s = 0;
            for k = 1:j-1
                s = s + L(i,k)*U(k,j);
            end
            L(i,j) = (X(i,j) - s)/U(j,j);
        elseif i < j
            s = 0;
            for k = 1:i-1
                s = s + L(i,k)*U(k,j);
            end
            U(i,j) = X(i,j) - s;
        end
    end
end
disp(L)
disp(U)
end